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Chapter 1 



1.1 Introduction 

As a result of uncontrollable interactions between quantum systems and their 
local environments, complex correlations develop between them which lead to 
the phenomena of decoherence and relaxation when only the quantum system is 
observed [l] [21 |3l |4] . As almost no quantum states can ever be completely iso- 
lated from their surroundings, the dynamics of so-called open quantum systems 
appear in almost all experiments in quantum physics, chemistry and biology, 
and a detailed understanding of the role of uncontrollable, noisy environmental 
interactions is required to extract genuine quantum effects from realistic data. 

In many cases, such as quantum optics and atomic physics, the effects of 
these processess are weak and relatively benign; while environmental interac- 
tions do degrade quantum effects, they do so on much slower timescales than 
those on which the effects operate and can be probed [5J [SI H] . Under these con- 
ditions, these quantum effects cannot just be unambiguously observed, they can 
even be controlled and potentially harnessed in new breeds of quantum device 
which can greatly outperform their classical analogues [J. Yet in other physical 
settings, such as the solid state and biological systems, the often strong and 
complex environmental interactions rapidly degrade quantum effects. Indeed, 
in many biological systems is has often been thought that relatively strong envi- 
ronmental noise is essential for directing an essentially classical and irreversible 
- c.f. reversible unitary dynamics- migration of energy through the complex 
energy landscape which connects energy producing and energy consuming parts 
of the system |l[i[10]. 

A good example of this latter paradigm is provided by pigment-protein com- 
plexes (PPCs) in photosynhetic organisms [51 [TU]. These structures are involved 
in the early stages of light-harvesting and excitation energy transport (EET) 
which initiate the carbon-fixing reactions of photosynhesis. The wide variety 
of PPC structures share the common motif that they contain optically active 
chromophore molecules embedded in a protein matrix which co-ordinates their 
spatial distribution. In typical photosynthetic organisms, the PPCs are arranged 
so that particular complexes (antennae) absorb photons via the creation of 
electronic excitations (excitons) on their chromophores, whilst other complexes 
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transport these excitations to reaction centers where electrons are released for 
photosynthetic chemistry [10]. The passage of excitations from generation to 
consumption is generally achieved through the existence of energy gradients in 
the potential landscape of the inter- and intra-complex chromophores [3 [9l [10] 
, which allows funneling of energy through dissipative processes induced by 
fluctuations of solvents and surrounding proteins. Remarkably, for many photo- 
synthetic systems under low light conditions, the quantum efficiency of photon 
capture, transport and charge generation is close to 100%P [9l [TO]. 

Although the migrating excitations in PPCs may be of a quantum mechan- 
ical nature, it was normally assumed that the complex, high temperature envi- 
ronments of functioning PPCs would rapidly destroy inter-exciton coherences. 
Consequently, the dissipative funneling of energy could be intuitively described 
and understood by effectively classical rate-equation dynamics such as those 
provided by the Forster and Dexter theories [H [9]. However, a much more 
complex picture of EET has recently emerged with the discovery of evidence 
for long-lasting inter-exciton coherences in the EET dynamics of the Fenna- 
Matthews-Olson (FMO) complex [11]. This complex is extracted from green 
sulphur bacteria, and functions like a biomolecular 'wire' that transports ex- 
citons from the light-harvesting chlorosomes to the charge-separating complex 
known as the reaction center [TOl [121 US]. Since the discovery of this evidence, 
similar effects have also been observed in complexes from marine algae and 
green plants [Ml [15] , and further FMO experiments have now suggested coher- 
ence lifetimes of around 1.5 ps at 77 K and a few hundred femtoseconds at 277 

K [iii[T7i[Ti[ini[5ni. 

These inter-exciton coherence times are striking, as they are almost an or- 
der of magnitude longer than the coherence times of single excitonic transitions 
(~ 100 — 200 fs) [17], and as a result they persist over a significant fraction 
of the total transport time in typical PPCs [l|. It has therefore been suggested 
that coherences may play an important role in driving the directed, highly- 
efficient EET observed in these complexes, and understanding this may provide 
valuable insights into how similar efficiencies could be achieved in artifical light- 
harvesting systems. However the mechanisms which preserve these coherences 
are currently unknown, and this and the intrincate interplay of noise and co- 
herence that generates efficient transport has become a very rich and active 
problem [22 [23l [24l [25l [26l [27l [Ml [22 ISQl [3ll [32 [33 

The dynamical behaviour of interacting open quantum systems is frequently 
investigated in terms of simple dynamical models in which environmental de- 
phasing and relaxation are treated with Lindblad or Bloch-Redfield master equa- 
tions. These methods are both based on the assumptions of weak system-bath 
coupling and the Markov approximation. However, these approximations are 
not valid for many realistic systems, and assuming that the correlation time 
of the environments in these systems is much faster than the system dynam- 
ics is frequently not justified. For instance, in typical PPCs the dynamical 



^ The transport time for a single excitation to pass throuh the FMO complex is estimated 
to be ~ 5 ps [S]. 
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timescales of the bath can be comparable or even slower than the EET dy- 
namics [35l |34j [30]. Moreover, in the limit of slow bath dynamics, pertur- 
bative treatments of the system-environment coupling cannot be used even if 
the system-bath coupling is intrinsically weak. Recently, important steps have 
been taken towards the development of non-perturbative and non-Markovian ap- 
proaches, including generalised approximate master equations [551155] . formally 
exact master equations which are unravelled by numerical hierarchy techniques 
(NHT) [35] |37], stochastic methods [38j, and numerical path integral (NPI) 
techniques such as quantum monte carlo [39] and QUAPI[30] [40]. There are, 
however, limitations concerning the quality of the uncontrolled approximations 
made in some approaches [33] [36] [38] , the restricted environmetnal structures 
accessible to several of these techniques [55] [57] , and many of the numericlly- 
exact methods are expected to become less efficient with decreasing tempera- 
tures [39] [30] [33 . Given the detailed information about the real protein spectral 
densities in PPCs is only just beginning to emerge [4T], a technique is required 
that can simulate EET for arbitrary spectral densities and coupling strengths, 
thus allowing experiments carried out under different conditions, including low 
temperatures, to be analyzed within one framework. 

The advanced numerical techniques such as NPI or NHT are approaches 
which deal with the time-evolution of the reduced density matrix of the quan- 
tum sub-system. The origins of the computational effort required to evaluate 
these schemes stems from the fact that without a separation of scales in PPC 
problems, the system and environment participate in the dynamics on an es- 
sentially equal footing. At a global level the dynamics has the character of a 
strongly-correlated many-body problem, suggesting an alternative approach to 
the problem based on numerical condensed matter theory methods. Because 
of the large number of environmental degrees of freedom, a direct simulation 
of the system and the environment appears rather daunting, but a number of 
powerful methods such numerical renormalisation group and sparse ploynomials 
space approaches have recently been developed to do precisely this [42] [43] [44] . 
The key to the success of these methods is that the dynamics of the system- 
environment space can be accurately reproduced in a truncated Hilbert space, 
which is intimately related to the fact that many standard open-system Hamil- 
tonians have an effectively ID structure which only contains nearest-neighbour 
interactions [45] . 

In this chapter we introduce another many-body appoach to open quantum 
systems simulation that combines an exact analytical mapping of the problem 
onto an effective ID nearest- neighbour model and the time- adaptive density 
matrix renormalisation group (t-DMRG) technique |46| . Since its introduction 
the t-DMRG technqiue has proven to be one of the most powerful, accurate 
and versatile methods for simulating many-body dynamics in 11? [46^, and in 
many cases leads to numerically exact results. The mapping we use to generate 
the ID representation also uses a novel application of the theory of orthogonal 
polynomials and a considerable portion of this chapter deals with this formalism 
and the physical interpretation that this alternative picture provides for open- 
system dynamics. 
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This chapter is organised as fohows: Section [1.2.11 introduces the standard 
open-system Hamiltonian and discusses the assumptions of this model. Sec- 
tion 11.2.21 sets out in detail the formal mapping technique that generates an 
equivalent ID representation of the open-system Hamiltonian which can be ef- 
ficiently simulated by t-DMRG methods. Section [1.2.31 points out a number of 
fundamental results on open-system structures which are revealed by this for- 
mal transformation, and points out how these might be used to increase the 
efficiency of future simulations. Section [1.31 presents numerical examples of our 
mapping and t-DMRG approach, one of which points towards a novel mecha- 
nism for long-lasting excitonic coherence in PPCs. Finally, a set of conclusions 
and future prospects for this approach are given in Section 11.41 

1.2 Open-system Hamiltonians and chain map- 
pings 

1.2.1 Standard model of open quantum system 

In this section wc shall consider the most common model of an open quantum 
system, in which a quantum sub-system interacts with a macroscopic number 
of environmental degrees of freedom and the total sub-system and environment 
state evolves under a purely unitary dynamics. Dissipation and decoherence 
appear when the sub-system is observed without any knowledge of the state 
of the environment, leading to a non- unitary effective dynamics for the sub- 
system's reduced density matrix. The total Hamiltonian can be written as H — 
Hs+Hj+Hb, where Hg is the Hamiltonian of the quantum sub-system's degrees 
of freedom, Hb is the free Hamiltonian of the environment and Hj describes 
the interaction of the system and bath variables. For the typical problems 
described in the introduction, the quantum subsystem consists in a finite number 
of quantum states which we denote as \i), and the system Hamiltonian can then 
be written in the general form 

N N 

i=i j=i 

where iJ is a Hermitian matrix and N is the total number of states which 
describe the quantum sub-system. For the excitation transport problems men- 
tioned in Section [1.1[ it is natural to associate the states \i) with the presence of 
an excitation on a physical, spatially localised site, in which case the diagonal 
matrix elements Ha give the local energies of these states and the off-diagonal 
matrix elements Hij quantify the probability amplitudes for these excitations 
to tunnel between sites i and j. 

For pigment-protein complexes and an extremely wide range of systems en- 
countered in physics, chemistry and biology, it is common to model the envi- 
ronment as a continuum of harmonic oscillators which interact linearly with 
the operators of the system [Ij jSj |42l [34] . We shall represent such an oscillator 
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environment in an explicit continuum representation [471 142] , which allows us 
to write Hj as 

N ^1 

Hi^y^Vs^ dkh,{k){a,{k) + alikj) (1.1) 

where Vsi are operators which act locally at site i and hi{k) describes the cou- 
pling to field modes labelled by a continuous quantum number k. The modes are 
described by creation and annihilation operators aj(fc) and ai{k) respectively, 
which obey the bosonic commutation relation, 

[a,{k),a]{k')]=S,jS{k-k'). 

We shall assume that k lies within the finite interval [0, 1], leading to an env- 
iornment with a sharp, finite bandwidth. Equation (jl.ip is not the most general 
form of linear system-bath interaction, for instance the bath(s) could couple to 
multiple operators at each site, couple to collective modes of the sub-system 
or have a non-linear interaction in the bath variables. The restricted form we 
use is motivated by the usual assumption in EET problems that the primary 
effect of the environment is to induce fluctuations of the local site energies Ha 
133] . We have also assumed in Eq. (jl.ll) that operators on each system site cou- 
ple to independent (commuting) environments, and we will therefore not deal 
with the issue of spatially-correlated fluctuations |1 7] |48] (361 l49l |4T| . The free 
Hamiltonian of the oscillators is 

Hb= f dkg,{k)a\{k)a^{k), (1.2) 
Jo 

where gi{k) is the dispersion of the field modes. The maximum frequency of the 
environment ujc is given by Wc = 5(1). The model is completed by specifying the 
spectral function of the environment J{i^), which in terms of the microscopic 
parameters of the Hamiltonian is given by [42], 

J(^u;) = h'[g-\c.)]^^^. (1.3) 

In Eq. (|1.3p . g^^{x) is the inverse function of the dispersion i.e. g^^{g{x)) = x. 
For all open-system problems where the environment is initially in a gaussian 
state, it can be shown rigourously that the influence of the environment on the 
reduced system dynamics is completely determined by J{lo) only [1] [3] l35j [2]. 
Equation (|1.3|) provides a relation for obtaining J{uj) from a specific miscroscopic 
interaction model, however it is often the case that the spectral function itself 
is given or assumed, in which case the functions h{k) and g{k) are not uniquely 
specified. In the following section we will work with a fixed J{uj) and use the 
indeterminacy of h{k) and g{k) to effect the mapping we shall now present. 
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(a) 




(b) 




J 



Figure 1.1: (a) Standard representation of a quantum dimer system in which 
each site is coupled to an independent continuum of harmonic oscillators, (b) 
After a unitary transformation of the oscillators, the entire system can be rep- 
resented as a ID chain with nearest neighbour interactions t„ and local energies 
e„. This equivalent many-body system can now be simulated efficiently using 
conventional t-DMRG techniques. 

1.2.2 Unitary transformation of the environment 

In this section we present the essential details of how we can convert the standard 
Hamiltonian structure of the open quantum system shown in Fig. Il.ll a into a 1- 
D form suitable for t-DMRG simulation. The key insight is that the interaction 
of the quantum sub-system with all the modes of environmental oscillators is 
equivalent to the local interaction of the sub-system with one end of an infinite 
ID chain of coupled harmonic oscillators as shown in Fig ll.ll b. The existence of 
a chain representation of the environment has been known for quite some time 
in a variety of quantum and classical contexts [J U^J [SOI [SIl ISH [S3], and has 
been of particular use in the study of quantum impurity problems by numerical 
renormalisation group methods ^42^ i43). In almost all previous approaches, the 
representation of the environment as a chain is used as an intermediate step that 
permits the application of a numerical technique. Consequently, the unitary 
transformation (see below) which maps the original open-system Hamiltonian 
onto a ID chain is often carried out numerically, following a discretisation of 
the continuous environmental spectrum to make the problem computationally 
tractable. However, these numerical mappings can often be numerical unstable, 
even for relatively unstructured environments. 

In our approach we carry out the mapping formally, using the theory of 
orthogonal polynomials to perform the mapping exactly and analytically. This 
formal approach allows us to make use of many of the rigourous results of 
orthogonal polynomial theory, and we shall show how their application reveals 
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universal properties of open quantum systems which are independent of the 
specific forms of the environmental spectral function. Orthogonal polynomials 
also have rigourous connections to other important mathematical objects, such 
as continued fractions, Cauchy transforms and random matrices, and our theory 
provides a very general framework for investigating how these objects might also 
be applied to the problem at hand. 

A vast literature on orthogonal polynomials exists, and research into orthog- 
onal systems is still extremely active, not least because of their important role 
in numerical quadrature, multi-dimensional interpolation, stochastic modelling, 
random matrices, approximation theory and analysis (SU [55l [56l [57l [58l [59] . 
A classic text is that of Szego [60], and many other fine books on the subject 
can be found in [SH IHIl [121 [Ml IM]- For the most part, the material presented 
in this chapter only deals with the simplest types of orthogonal polynomial on 
the real line, and in what follows we shall use several standard results without 
proofs. The detailed proofs can be found in the any of the books above, but are 
also conveniently collected together in the context of the open quantum system 
problem in Chin, Rivas, Huelga and Plenio |65) . 

The starting point of the mapping is a unitary transformation which acts 
just on the environment oscillators. In order to prevent too many subscripts 
and summations from cluttering up our presentation, we shall only consider a 
single system site in what follows, dropping the site index i throughout. As 
our open-system model consists of independent baths coupled to each site, the 
extension to multiple sites is trivial, and will be touched on again briefly in 
Section 11.41 We implement the transformation by defining new bosonic modes 
according to 



where h(k) is the coupling function in Ec^. (|1.1|) . 7r„(fc) is a nth monic or- 
thogonal polynomial (to be defined below), and p„ is a normalisation constant. 
The corresponding transformation for 6jj is obtained by taking the Hermitian 
conjugate of Eq. (|1.4p . and we note here that the parameters /i(fc),7r„(fc) and 
Pn are all real- valued. The function 7r„(fc) is a monic nth degree polynomial 
T^nik) = ^^^oCnjk^ where the monic condition means that c„„ = 1. The 
coefficients of the polynomials cj are chosen so that they obey the following 
orthogonality condition, 



which defines the normalization constant appearing in Eq. ()1.4|) . The poly- 
nomials 7r„ are known as monic orthogonal polynomials (MOPs) of the weight 
function /i^(fc). For a strictly positive weight function, as is manifestly the case 
for the weight fmiction h'^{k), a complete sequence of MOPs can always be 
found as a result of Favard's Theorem [601 [£l]- The orthogonality condition 




(1.4) 




(1.5) 
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immediately implies that, 

[bn,blj = pnPm dk dk'h{k)h{k')'Kn{k)Trm(k')[a{k),a\k')] 



= PnPni / dkh {k)T:„{k)'Km{k) 

Jo 

= 5nm, (1-6) 

where we have used the commutation relation of the continuum field modes 
in the second line and the orthogoanlity relation of Eq. ()1.5|) in the third. 
The transformation is real orthogonal and preserves the bosonic commuation 
relations of the new modes bj^ H. The inverse transformation is given by, 

oo 

= X! Hk)Pn'^n{k)bn, (1.7) 

which we now use to construct the chain Hamiltonian by substituting Eq. ()1.7|) 
into the original open-system Hamiltonian H — Hs + Hj + Hb ■ The transforma- 
tion of the environment modes does not affect the system operators, and there- 
fore Hs and the system operator Vs in the interaction term Hj are unchanged 
by this operation. Let us now consider the effects of the transformation on 
the interaction Hamitonian Hj and free bath Hamiltonian Hb separately. The 
interaction term Hj transforms in the following way. 

Hi = Vs f dkh{k){a(k) + a\k)) 
Jo 

= VsYpnibn + bi) dke{k)nn{k) 



n=0 ""^ 
oo „i 

= VsYpnibn+bi) dkh^{k)7ro{k)TTn{k) 

n=0 

= Vsp^\b^ + bl), (1.8) 

where we have used the fact that - by definition - 7ro(fc) = 1, and then the 
orthogonality relation in the last line. The result of the transformation is that 
the system now couples to only a single mode 60 of the new representation of 
the environment. We now turn to the bath Hamiltonian Hb- This transforms 
into, 

Hb = [ dkg{k)a\k)a{k) 

dkg{k)TTn{k)TT„^{k). (1.9) 



00 00 ^ 1 

E E blMn / 
n=Om=0 "^0 



^ We note that this transformation, and everything that follows in this section, would also 
hold true for an environment of fermionic oscillators. 
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At this point we cannot proceed further until the dispersion function g(k) is 
specified. As discussed in Section 11.2.11 the open-system dynamics of the sub- 
system is completely determined by the spectral function. Therefore, for a given 
J{(jj) we have the freedom to choose the form of either h(k) or g{k) as long as 
Eq. (jl.Sp is respected. For reasons that will soon become apparent, we choose 
to take the dispersion to be g{k) — Wcfcd. This fixes /i^(fc) = LUcJ{uJck), and thus 
the MOPs of our transformation are orthogonal w.r.t. a weight function which is 
just proportional to the spectral function. The freedom to partition the spectral 
function between h(k) and g{k) is also used in the NRG approach [42]. where 
it is used to logarithmically discretize the spectral function. The subsequent 
mapping onto a chain can also be solved analytically with generalised MOPs, 
and an example of such a solution is given in section 11.2.51 The linear form of 
g{x) now allows us to use another general property of MOPs, which is that they 
all obey the following three-term recurrence relation |60j |65] , 

kTTn{k) = a„7r„(/c) + ^„7r„_i(fc) + 7r„+i(A;), 7r_i(A;) = 0, (1-10) 

where the sequence of numbers a„ , /3„ are unique for a given weight function and 

are given by a„ = dfc/i^(fc) fc7r„(fc)7r„(fc) and /3„ = p„p„+i dfc /i^(fc)fc7r„(fc)7r„_i(fc). 

If we now substitute the linear form of ^(fc) into Eq. (|1.9p and use the recurrence 

and orthogonality relations, we obtain 




n— m— 





(1.11) 



Due to the choice of the linear dispersion, the transformed bath Hamiltonian 
takes the form of a one dimensional harmonic chain with only nearest neighbour 
coupling. From the definitions of /3„ and p„ one can easily show that Pn/ Pn+i = 
Pn+i, allowing us to rewrite Hb in the final, symmetrised form, 



oo 




(1.12) 



n=0 



where e„ = ujc<^n and t„ = u}c\/ Pn+i- We have now completed the formally ex- 
act transformation from the original Hamiltonian to a 1 — I? nearest-neighbour 



•^One can make other choices for the dispersion and this actually allows a number of different 
types of chain structures to be generated. The potential uses of these generalised structures 
in numerical applications is an interesting and open topic. 
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Hamiltonian. Collecting together all the transformed terms, the total Hamilto- 
nian in the chain representation is given by, 



Htotai =Hs+ rjVsibo + bl) + Y, ^nb^hn + t^hlhn+i + ^„foi+i&„, (1.13) 
where using h'^{k) — ujcJ{uJck) we have defined the coupling constant ry, 

dkh^{k) = / dujJ{uj). (1.14) 
Jo 

The dynamics of the many-body system and bath state under this Hamilto- 
nian structure can now simulated using t-DMRG, as will be described in Section 
11.31 However before presenting the simulation technique we shall briefiy describe 
some physical implications of the exact mapping. 

1.2.3 Universal properties of continuous environments and 
the determination of the chain frequencies and cou- 
phngs 

In section 11.2.21 we derived the relation between the chain oscillator frequencies 
e„ and the couplings to the MOPs recurrence coefficients a„ and /?„. The 
determination of the chain which corresponds to an environment characterised 
by a given J(a;) therefore reduces to the problem of determining the recurrence 
coefficients of the MOPs w.r.t. the weight function J{uj). For several important 
weight functions, these recurrence coefficients can be given by simple analytical 
formula. A comprehensive list and analysis of these classical MOPs can be found 
in [HJ EOl EU |62l [63l 111] . A very useful example is provided by the shifted- 
Jacobi polynomials P°'^(fc), which are defined on the interval k G [0, 1]. These 
polynomials are orthogonal w.r.t. the Caldeira-Leggett spectral density /weight 
function J[uJck) = aoj].~''uj'^9{ijJc — ^^), which is often used in discussions of the 
spin-boson model, the archetypical model of an open quantum system [1] [31 [SSI 
[3H [331 mi [SZ] ■ The corresponding chain frequencies, inter-chain couplings and 
coupling to the quantum system are given by. 




UJr 



2 V (s + 2n)(2 + s + 2n) 



(1.15) 



ujc{l + n){l + s + n) ji + s + 2n 



(s -t- 2 + 2n)(3 + s + 2n) V 1 + s + 2n 

- (1.17) 

1 + s 

As one can see from Eqs. (Il.isp . (|1.16l) and (|1.17p . the energy scale of the 
bath Hamiltonian and interaction terms are set by Wc (as one would expect) 
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(b) i Ky< >*KV --- 



Figure 1.2: Illustrative sketch of open-system dynamics in the chain represen- 
tation, (a) Sub-system initially injects excitations (shown as wavepackets) into 
inhomogeneous region of the chain. Scattering from inhomogeneity causes back 
action of excitations on the system at later times and leads to memory effects 
and non-Markovian sub-system dynamics, (b) At long times, after multiple 
scattering, excitations penetrate into the homogeneous region and propagate 
away from the system without backscattering. This leads to irreversible and 
Markovian excitation absorption by the environment. 



and the total dynamics of the open quantum system are determined by the 
dimensional-less parameters a and the eigenvalues of Hg/ujc- In addition to 
this, we can also immediately infer the asymptotic {n — > oo) parameters of the 
chain, finding that e„ — > Wc/2 and t„ — )■ as n — > cx). These asymptotic 

values do not depend on the values of s which characterise the shape of J{lo) and 
are thus universal for spectral densities of the Jacobi power-law form. At large 
distances from the sub-system the harmonic chain becomes homogenous and 
excitations in this part of the chain become simple harmonic waves. Using the 
asymptotic values of e„ and tn one can simply diagonalise the homogenous part 
of the chain, yielding the dispersion VL{q) — \^c{^ ~ cos(7rg)) of excitations with 
wavevector q. As shown in Fig. 11.21 the asymptotic region of the chain can be 
loosely thought of as a type of 'transmission' line, whose homogeneity ensures no 
backscattering of excitations towards the system. As sketched in Fig. (|1.2p . this 
enables this region to carry away excitations from the sub-system irreversibily 
at long times, as one would expect for a dissipative environment. On physical 
grounds we would also require that this region of the chain should be able to 
support excitations at all frequencies covered by the original spectral function, 
and indeed it can be seen that the asymptotic values e„ and i„ are the only 
values which give the correct bandwidth for the asymptotic region. 

The emergence of a universal asymptotic chain appears directly from the 
analytical formula for the Jacobi recurrence coefficients, and can be physically 
motived by the arguments given above. Indeed, on the basis of the physical 
arguments one might expect this asymptotic homogeneity to appear for any 
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finite bandwidth environment, and this indeed turns out to be the case. The 
proof is due to Szegcjl, who was able to show that the asymptotic values for 
Cn and tn are, respectively, Wc/2 and Wc/4 for any weight function h'^{k) which 
obeys the inequality [601 ES] 



Weight functions which obey Eq. (|1.18p are said to belong to the Szego 
class [61] [54| . In the context of the open-system problems we are considering, 
a huge range of spectral functions over a finite interval falls within the Szego 
class and thus the homogenous asymptotic chain appears in almost every chain 
representation of a physical, finite bandwidth environment. Notable example of 
non-Szego spectral function corresponds to spectral functions containing band 
gaps or spectral functions defined over semi-infinite domains. We shall not 
consider these cases in this chapter, but they are dealt with in Chin, Rivas, 
Huelga and Plenio [65]. 

The existence of a universal asymptotic form of the chain region leads to 
a very appealing and simple picture of memory effects and non-Markovianity 
in open-system dynamics. The chain structure itself implies a natural causal- 
ity, or set of timescales, over which different regions of the chain contribute 
to the dynamics as shown in Fig. (|1.2|) . At early times the system interacts 
with the modes on the left of the chain, injecting excitations into this region 
which then begin to propagate to the right. Because of the inhomogeniety of 
this region, which is dependent on the specific form of the spectral function, 
these excitations will undergo scattering and some of them will return and act 
on the system at a later time. These backscattering proceses represent memory 
effects in the system-bath interaction and depend sensitively on the form of the 
spectral function. At later times excitations propagate into the homogeneous 
asymptotic region of the chain and are effectively absorbed irreversibly by the 
environment. The dynamics of this process is independent of the shape of the 
spectral function and describe a dissipative, long-time Markovian dynamics of 
the sub-system. Therefore in the chain representation the bath corrrelation 
time and related memory effects are associated with the typical time it takes 
an excitation to exit the inhomogenous region close to the system. This time 
depends on the form of the bath which determines the size and spatial ex- 
tent of the backscattering potential seen by these excitations. The strengh of 
non-Markovian effects on the sub-system dynamics depends on how excited the 
inhomogenous region is during the time evolution and will therefore depend on 
the rate at which excitations are injected into this region i.e. it will be depen- 
dent on the coupling strength. Non-trivial, initially non-Markovian dynamics 
is therefore expected when the excitation injection rate is much larger than the 
rate of escape from the inhomogenous region of the chain at early times. 

*We have rephrased the theorem in terms of the chain parameters. The theorem pre- 
sented in Refs. [601 165) is actually a statment about the asymptotic values of the recurrence 
coefficients of a sequence of orthogonal polynomials defined over a finite interval. 




(1.18) 
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(a) 



(b) 




Figure 1.3: (a) Diagonalising the homogeneous part of the chain after site N 
leads to an effective environment acting on this terminal oscillator, as shown in 
(b). This terminal spectral density is universal for any spectral density in the 
Szego class, suggesting that complex environments may be efficiently handled 
by only treating the initial oscillators of the chain which encode the specific 
characteristics of a given environment. 
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Another representation of this idea of a non-trivial, non-universal early time 
dynamics which evolves into a universal dissipative dynamics is shown in Fig. 
(jl.3p . After the chain parameters have, to within some arbitrary tolerance, con- 
verged to the asymptotic values at site N, the remainder of the chain to the right 
is diagonalised. This provides an effective environment acting on the A'^th mem- 
ber of the chain which in the limit N oo possesses a universal spectral function 
proportional to the Wigner semi-circle distribution Jt{'^) oc ^Juj{uJc — w), which 
is an important equilibrium distribution in random matrix theory and which is 
also imitately related to the properties of Chebyshev polynomials [65l [571 HI] • 
This representation suggests a possible reduction in the complexity of simulat- 
ing the dynamics of a complex environment, as in many cases the convergence of 
the chain parameters is rather rapicd. It may therefore be possible to simulate 
the system by treating only the first few non-trivial sites of the chain explicitly, 
and then using numerically cheaper semi-classical, or even classical, methods to 
model the damped mode at site A^. 

As t-DMRG simulates the entire wavefunction of the system and environ- 
ment, we will be able to explore the correlations and entanglement between the 
system and bath, allowing us to accurately assess the quality of such an approx- 
imation and how to improve upon it systematically. Investigating system-bath 
correlations may also be of some relevance for understanding how entangle- 
ment is generated between different components of open systems [26l [27l |32] , 
and is of direct relevance for the recently-developed theory of measures of non- 
markovanity [68l [69] . An important practical application of having access to 
bath information is that we can also explore at the microscopic level how prepa- 
ration and propagation of wave packet dynamics in complex environments can 
influence EET networks. 

This idea of the reduction of complex environmental spectra has also been 
addressed by Burghardt et al. [HI [52l [53] , who have derived an iterative formula 
for the effective spectral density acting on site A^ as A^ is increased. Using a 
mass-weighted co-ordinate representation of the environment and chain, Mar- 
tinazzo et al. also empirically found that the spectral density converges to a 
universal limit under certain conditions, and that this terminal spectral den- 



sity has the Ohmic Rubin model form Jt{'^) oc lo Jl - ^ [iS]- Their method 
makes extensive use of continued fractions and Cauchy tranforms, which are in- 
timately related, via the Jacobi matrix, to orthogonal polynomials [54 ] l60 l l6T|. 
As shown in Weiss [T], the Rubin spectral density can also be represented by a 
coupling to a uniform chain of harmonic oscillators coupled by nearest neigh- 
bour interactions, and the formal links between these approaches is currently 
being investigated within the framework of orthogonal polynomial theory. 



^For the Jacobi spectral functions the parameters converge to their asymptotic values as 
^ as n — > oo. 
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1.2.4 Continuous, discrete and mixed spectral densities 

In the previous section we dealt with chain representations related to continuous 
spectral functions over a finite interval. In many situations we also encounter 
spectral densities containing discrete contributions, either as a result of the 
physical presence of strong coupling to discrete modes of the environment or an 
artifical discretisation of the environment that has been performed to facilitate a 
numerical approach to the problem. As discussed in Refs. [54l [60l |6T1 l63l l62l [65] 
it turns out that MOPs can also be found for such spectral functions, permitting 
the formal transformation of these problems into the ID harmonic chain prob- 
lem. We shall illustrate this with analytical results for the important case of 
a logarithmically-discretized power-law spectral density. This artifically discre- 
tised spectral density plays an important role in the powerful numerical renor- 
malisation group approaches to quantum impurity problems [42]. Numerical 
results for a physical spectral density with a discrete component will be pre- 
sented in Section fl. 3. 2 1 

To handle discrete components we consider spectral functions J(fc) of the 
form, 

JV 

J{k) = h^{k) + ^wjS{k- kj) k,kje [0,1], (1.19) 

where h^{k) is a continuous, non- negative spectral density, Wj are positive 
weights for discrete contributions to the sepctral density and kj the (scaled) 
frequencies at which these discrete features occur. Under these conditions it 
can be shown that a set of MOPs can always be found which obey [54] [70] , 

1 nl N 

dkJ{k)TTn{k)Trjn{k) = / dkh'^{k)Trn{k)n„i{k) + y2wj7rn{kj)Trm{kj) 
Jo ^.^1 

= Pn^Snm, (1.20) 

and that these MOPs possess the key properties we need to implement the 
chain transformation, such as the three-term recurrence relation. In the extreme 
case where all wj's are zero, an infinite sequence of MOPs, like those we have 
already considered, arises. In the opposite extreme where h?{k) — 0, there 
is a finite number N of discrete MOPs which obey the discrete orthogonality 
condition X^jLi Wjnn{kj)TTm{kj) = Snm- Just like in the continuous case, there 
exists a number of classical discrete MOPs whose properties can be expressed 
in analytical form, and a comprehensive list can be found in [63j . In the mixed 
case the sequence of MOPs is also infinite, and while a few special cases can be 
solved analytically [54] [61] , the MOPs for these cases normally have to be found 
numerically. 

For the general mixed spectral density of Eq. ()1.19|) a number of very ef- 
ficient algorithms have been developed for computing the values of the recur- 
rence coefficients a„,/3„ which enter the chain Hamiltonian. The most effective 
of these for mixed problems involve adaptable discretisation and quadrature 
schemes which are collected in W. Gautschi's software package ORTHOPOL 
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[701 . These algorithms were used to determine the chain parameters for the 
numerical t-DMRG results in Section [T73l 

Before presenting numerical simulations we shall quickly give a practically 
useful example of an analytical solution to an important and purely discrete 
MOP problem. 



1.2.5 Logarithmically-discretised spectral density 

In numerical renormalisation group (NRG) studies of quantum impurity prob- 
lems of the spin-boson model- type [47l |42l |43] , a Hamiltonian of the form of 
Eq. (|1.2.1|) is first logarithmically discretized following the procedure set out 
in [42] . The Hamiltonian H after the logarithmic discretisation of the reservoir 
continuum takes the discrete form 

^ oo oo 

H = Hs + —j= ^ 7„(a„ + at ) + ^ C«aUn: 

"'V I" „^o n=0 

where 



1 + s ^ ^ ^ ^ 

Cn = — ^ — — a;,A-". (1.22) 

It has been shown by Bulla et al. that this Hamiltonian can then be mapped 
to a nearest- neighbour chain Hamiltonian of the form [4 7) . 



H, = H, + ly'^Vsibo + bl) + ^w„6t 6„ + t„6l+i6„ + i„6j,6„+i, (1.23) 

n 

by a real orthogonal tranformation bl^ = Unm<An provided that the matrix 
elements Unm obey the three-term recurrence relation, 

CnUmn — ^raUran "t" tvaUjri+ln ^ ^'m— l^m— In- (-^■^^) 

In the NRG approach, this recurrence relation is solved numerically by a sim- 
ple iterative procedure that rapidly becomes unstable as the size of the chain 
increases. However, the appearance of a real symmetric three-term recur- 
rence relation suggests that a closed form solution exists in terms of suitably 
chosen orthogonal polynomials. The resulting polynomials are in fact well 
characterised, allowing us to find the chain parameters of the logarithmically- 
discretised chain exactly. These polynomials are the little-q Jacobi polynomials 
p„(A^™, A"**, 1| A~^). These are normally not part of the classical scheme 
of discrete orthogonal polynomials and are in fact q-analogues of the classical 
Jacobi polynomials. A detailed discussion and a list of the other known q- 
orthogonal polynomials can be found in [71] . Their important properties for our 
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purposes is that that they obey the orthogonaUty relation, 

oo 

Sn,nN?, = ^A-^(l+^)p„(A-^A-^l|A-l)p„(A-^A-^l|A-lx,l•25) 

and the recurrence relation 

A->,(A-^A-^l|A-l) = (A,+C,)p,(A-^A-^l|A-l)-A,p,+l(A-^A-^l|A-l) 

- Qp,_l(A-^A-^l|A-l). (1.26) 

The normalisation constants Nn in Eq. (I1.25P and the recurrence constants 
An , Bn and C„ can be expressed in closed form and can be evaluated easily 
without any need for potentially unstable iterative techniques. The various 
coefhcients are listed in [65l [71]. As shown in Ref. [65], with just these two 
properties one can prove that the unitary matrix, 

_ A-"^Pn(A"",A-M|A-i) 

solves the recurrence relation of Bulla, and thus carries out the mapping exactly. 
The resulting chain parameters of Eq. (|1.23p are then given by. 



Cs{An+Cn), (1.28) 

Nn 



<s{^-]An. (1.29) 



1.3 Numerical results and applications 

We now demonstrate the implementation of our joint mapping and t-DMRG 
approach with some specific spectral densities of relevance for PPCs in photo- 
synthetic organisms. Although the richly structured environments used in the 
PPG literature are taken as challenging examples, it should be emphasized that 
this new simulation tool is completely general, and can be applied to any sys- 
tem linearly coupled to bosonic or fermionic environments of arbitrary spectral 
density. The PPG results were first presented in Prior et al. [72]. We shall 
consider a dimer system consisting of two sites 1 and 2 with local site energies 
ei, £2, and which are connected by a tunneling amplitude J. Each site is coupled 
to its own independent environment as in Eq. (jl.ip . and each environment is 
described by identical spectral densities. After the chain transformation, the 
system structure is exactly as shown in Fig. (|l.ip b. The initial state of the 
system for all simulations is taken as the separable pure state p = ps ® PB, 
where ps describes an initial excitation on site 1 and /Os is the vacuum state 
for the chain. The pure state initial condition implies that we are considering 
the open-system dynamics at zero temperature. These conditions on the spec- 
tral densities and states were chosen for simplicity and for their correspondence 
to the physical conditions found in PPGs immediately after photoexcitation. 
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but these conditions are not required for the successful implementation of our 
method. 

The pure state t-DMRG algorithm employed is the standard one presented 
in Refs. [ISl |731 [TH [7S] , which is used to evolve the total wavefunction of 
the dimer and environments in real time. Observables of the sub-system Os at 
time t were obtained from the expectation values {Os){t) = (^'(t)|Os|^(i)). In 
all t-DMRG simulations, we found that the results converged to less than 0.1% 
with just 11 bosonic levels per site, 30 Schmidt coefficients and 100 chain sites 
over the whole dynamics [7^ . 



1.3.1 The overdamped Brownian oscillator spectral den- 
sity 

To start with, we look at the overdamped Brownian oscillator spectral density 
which has been extensively studied in the context of PPG dynamics, and which 
can be solved in an numerically exact way in the high-temperature limit [35j . In 
our notation, the overdamped Brownian oscillator spectral density has a simple 
Ohmic form, 

= (1.30) 

where A is the reorganisation energy of the bath, defined hy X — J^" J{Lu)uj~-^dLu, 
and is taken as our measure of the site-environment coupling strength. The pa- 
rameter 7 approximately sets the dynamical response time of the bath, and the 
following simulations use values of 7 which are smaller than the dimer energy 
scales in order to observe non-Markovian effects 25> l30j i38j . 

Figure (|1.4[) shows the population on site 1 as a function of time for various 
values of A. For A < lOOcm"-'^ we find damped oscillations which persist for at 
least 1 ps. For larger A, coherent dynamics are always seen for a few hundred 
femtoseconds before the dynamics becomes incoherent, although as A increases 
the duration of coherent motion becomes shorter. For A > 200cm~^ the incoher- 
ent relaxation rate decreases dramatically, and an increasingly large population 
is trapped on site 1 over the timescale of the simulations. This quantum-Zeno- 
like phenomenon may be related to the well-studied localisation transition found 
in Ohmic and sub-Ohmic spin-boson models atT^OK [H l3l [76 l [77 l l44l l66] . 
This is a non-perturbative feature of the dynamics, and similar dynamics have 
also recently been observed in NRG and NPI studies of the sub-Ohmic spin- 
boson model l43l l67l. 



1.3.2 Other spectral densities 

We now demonstrate the versatility of our method w.r.t. the microscopic 
system-bath interactions by considering a much more complex and structured 
environmental spectral function taken from a recent study of photosynthetic 
EET. In Ref [21], Adolphs and Renger use a combination of super-Ohmic den- 
sities and a coupling to a single effective high-energy mode to model the envi- 
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Figure 1.4: Evolutions of the population on site 1 for the spectral density of Eq. 
(I1.30p at r = K, and various reorganisation energies A. Simulation parameters 
are J — lOOcm^^, ei — 62 = lOOcni^^ and 7 = 53cm^^ 
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ronment. In our notation this spectral function can be written as, 



2nX [lOOOw^e ( -"i ) ' + A.Suj^e ( 
9!(1000wf + 4.3a;f) 



) 



(1.31) 



where we have kept the relative contributions of the two continuous parts of the 
spectral density as they are in [21J , but have also introduced an overall reorgani- 
sation energy A to be used as a free parameter. The coupling to the high-energy 
mode is fixed, and the parameters of the simulation are J = lOOcm"^, ei — 62 = 
100cm~^,a;i = 0.5cm~^,a;2 = 1.95cm~^, w// = 180cm~^,uJc = lOOOcm^^ and 
Sh = 0.22 [51]. With these values the continuous part of J{uj) extends over a 
frequency range of about 900cm~"'^, and ujh is almost resonant with the energy 
difference ( 224cm~'^) of the dimer eigenstates oi Hg as the coupling strength of 
this mode to a site is 84cm^^. 

The chain transformation and DMRG method offers numerical advantages 
over some other techniques for spectral functions which contain delta functions 
or damped resonances, as strong coupling to such modes of the environment do 
not have to be considered as part of the system Hamiltonian. As we discussed 
in Section ll.2.4| discontinuous features in the spectral density simply modify 
the MOPs of the chain mapping, allowing simulation of an arbitrary number 
of such discrete mode interactions in the presence of a continuous background 
without any increase in the complexity of the simulation. Coupling to undamped 
modes with frequencies comparable to or smaller than the dimer energies have 
to be considered as part of the system in approaches like NPI, or if included in 
the spectral function, they must be artifically damped so that their long-time 
correlation function decays fast enough to be treated accurately within the finite 
memory time imposed on these methods. 

The interaction with the near-resonant oscillator has a pronounced effect 
on the population dynamics, and Fig ll.51 shows how this coupling leads to a 
coherent beating effect which periodically suppresses population oscillations for 
A < SOOcm^"'^. These coherent multi- frequency effects are a strong sign that even 
though we treated the discrete mode as part of the environment, the mapping 
and t-DMRG method accurately treat the quantum coherent interations with 
this mode. In situations where site 2 might transfer population to another 
system, such a coherent suppression of oscillations could lead to an enhancement 
of EET from the dimer to that system. As A increases, the continuous part of the 
spectral density dominates the dynamics and we observe qualitatively similar 
behaviour to the dynamics obtained in Fig. 11.41 We note that the trapping- 
like dynamics for large A is less severe for this super-Ohmic J{uj), although the 
dynamics are still highly non-Markovian for strong coupling. 

A particularly striking feature of Fig. 11.51 is that in the regime of optimal 
EET (A ^ lOOcm"^), the high-energy mode leads to low amplitude oscillations 
which persist for at least 1.5 ps. When the high-energy mode is decoupled, 
coherent oscillations vanish for A = lOOcni^^ after just 0.3 ps. Experimen- 
tal observation of such persistent undamped oscillations after a fast population 
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Figure 1.5: Evolutions of the population on site 1 for the spectral function of 
Eq. (|1.3ip at various reorganization energies A and T = K. Dimer parameter 
are J = 100cm~^,ei — £2 = lOOcni"^. Dashed line shows how the dynamics 
when the high-energy mode is decoupled. 

transfer could thus indicate the presence of discrete high-energy modes in the 
environment of PPCs, and could be a useful signature for determining realistic 
J(cj)s in these complexes [721 l2H]- We also note that broadening the discrete 
mode by replacing the delta function in Eq. (jl.31l) with an appropriate line- 
shape for a damped oscillator leads to damping of these long-lasting oscillations 
(not shown), indicating that these features are induced by the quantum nature 
of the interaction to the resonant discrete mode. Recent experiments on the 
FMO complex have observed extremely long electronic coherence times of 1 — 2 
ps, which could be consistent with the effects described above, as vibrational 
coherences are typically much longer-lasting than electronic coherences. 

1.4 Conclusions and future developments and 
applications 

In this chapter we have presented the formal development of a mapping tech- 
nique which converts the standard representation of open-system Hamiltonians 
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Figure 1.6: (a) A multi-site configuration with independent baths which could 
in principle be simulated using recent developments in t-DMRG techniques, (b) 
Multiple sites coupled in a correlated way to a common environment require the 
treatment of long-range interactions between the sub-system and chain. 

into a ID chain Hamiltonian with nearest-neighbour interactions. Using or- 
thogonal poynomials we have found a way to carry out this transformation 
exactly, and in doing so have rigourously demonstrated a number of hitherto 
unrecognised universal properties of typical open-system structures. Although 
this chain mapping is a fascinating subject in its own right, and one that is cur- 
rently being actively investigated, it also provides a representation which allows 
the powerful t-DMRG algorithm to be used in simulating open-system dynam- 
ics under complex, non-perturbative and structured environmental interactions. 
Such environments are thought to play an important role in photosynthetic ex- 
citation dynamics and the accuracy and versatility of the t-DMRG approach 
has been illustrated in our numerical examples, where it was discovered that 
discrete resonances in the spectral function can induce long-lasting coherent dy- 
namics of similar duration to those observed in some PPG complexes. Because 
this approach simulates the entire many-body wavefunction, it should also allow 
us to study the dynamical generation of correlations and entanglement between 
the system and bath, permiting us to explore the ideas of universality and bath 
reduction schemes presented in Section 11.2.31 An important practical applica- 
tion of this bath analysis would also be to examine in microscopic detail how 
vibrational wavepacket dynamics generated by sudden photoexcitation can ef- 
fect EET dynamics. The microscopic nature of the quantum states leading to 
the long-lasting electro-vibronic coherences can also be inferred from such an 
analysis. 

However, detailed simulations of the PPG systems which could be compared 
to experimental data require a number of technical developments of the method 
used to produce the results of Prior et al. [72] . The most obvious is the need 
to account for finite temperatures, a problem which have already been resolved 
with the recent development of mixed-state t-DMRG algorithms 78] . Another 
development is the extension of the method to multi-site networks with inde- 
pendent enviroments. Performing the chain transformation on such a system 
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leads to the Hamitonian structure shown in Fig. (I1.6l) a. This system can still 
be treated as as an effectively ID chain with larger local dimensions, allow- 
ing standard t-DMRG to be applied. Finally, many current theories about 
the long-lasting coherence in PPCs invoke the idea that spatial correlations 
of environmental fluctuations may support long-lasting quantum coherences in 
these structures. Assuming that the sites couple in different ways to a com- 
mon environment, the effects of spatial correlations can be investigated. The 
chain representation of such an open-system is shown in Fig. (jl.6p b. Although 
the number of environmental degrees of freedom to simulate is reduced, one 
is now faced with having to deal with longer-range interactions between the 
sites and chain. Extensions of t-DMRG to handle long-range interactions have 
also recently been developed. Taken together, these developments indicate that 
extremely efficient, accurate, and completely general algorithms for simulating 
open-system dynamics have come a step closer. 
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